Particle sizes of PVC plastic produced by a machine are measured. The machine is operated by three different people, and eight different batches of resin are used. Two measurements are made for each combination of these two experimental factors.
Main effects
pvcfit1 <- lm(psize ~ operator + resin, data=pvc)
summary(pvcfit1)
##
## Call:
## lm(formula = psize ~ operator + resin, data = pvc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9500 -0.6125 -0.0167 0.4063 3.6833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.2396 0.5226 69.345 < 2e-16 ***
## operatorBob -0.2625 0.4048 -0.648 0.52059
## operatorCarl -1.5063 0.4048 -3.721 0.00064 ***
## resinR2 -1.0333 0.6610 -1.563 0.12630
## resinR3 -5.8000 0.6610 -8.774 1.13e-10 ***
## resinR4 -6.1833 0.6610 -9.354 2.11e-11 ***
## resinR5 -4.8000 0.6610 -7.261 1.09e-08 ***
## resinR6 -5.4500 0.6610 -8.245 5.46e-10 ***
## resinR7 -2.9167 0.6610 -4.412 8.16e-05 ***
## resinR8 -0.1833 0.6610 -0.277 0.78302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.145 on 38 degrees of freedom
## Multiple R-squared: 0.8595, Adjusted R-squared: 0.8262
## F-statistic: 25.82 on 9 and 38 DF, p-value: 1.474e-13
confint(pvcfit1)
## 2.5 % 97.5 %
## (Intercept) 35.181635 37.2975319
## operatorBob -1.081983 0.5569834
## operatorCarl -2.325733 -0.6867666
## resinR2 -2.371544 0.3048775
## resinR3 -7.138211 -4.4617892
## resinR4 -7.521544 -4.8451225
## resinR5 -6.138211 -3.4617892
## resinR6 -6.788211 -4.1117892
## resinR7 -4.254877 -1.5784559
## resinR8 -1.521544 1.1548775
This model assumes the influence of the two factors is additive, the model only contains the “main effects”. The meanings of the coefficients are:
- “(Intercept)” is particle size for Alice and R1
- “operatorBob” is particle size for Bob relative to Alice
- “operatorCarl” is particle size for Carl relative to Alice
- “resinR2” is particle size for R2 relative to R1
- “resinR3” is particle size for R3 relative to R1
- (etc)
We can use anova( ) to test if there is evidence either of these main effects is important. For example, to test if there is evidence that the operator is important to the outcome we can test pvcfit1 against a model in which operator is dropped:
pvcfit0 <- lm(psize ~ resin, data=pvc)
anova(pvcfit0, pvcfit1)
## Analysis of Variance Table
##
## Model 1: psize ~ resin
## Model 2: psize ~ operator + resin
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 40 70.533
## 2 38 49.815 2 20.718 7.902 0.00135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Contrasts and confidence intervals
anova( ) lets us test if a particular factor or interaction is needed at all, and summary( ) allows us to see if any levels of a factor differ from the first level. However we may wish to perform different comparisons of the levels of a factor – this is called a “contrast”. We might also want to test some more complicated combination of coefficients such as a difference between two hypothetical individuals. In general this is called a “linear hypotheses” or a “general linear hypothesis”.
Say we want to compare Bob and Carl’s particle sizes. We will use the pvcfit1 model.
coef(pvcfit1)
## (Intercept) operatorBob operatorCarl resinR2 resinR3 resinR4
## 36.2395833 -0.2625000 -1.5062500 -1.0333333 -5.8000000 -6.1833333
## resinR5 resinR6 resinR7 resinR8
## -4.8000000 -5.4500000 -2.9166667 -0.1833333
K <- rbind(Carl_vs_Bob = c(0, -1,1, 0,0,0,0,0,0,0))
K %*% coef(pvcfit1)
## [,1]
## Carl_vs_Bob -1.24375
This is the estimated difference in particle size between Carl and Bob, but can we trust it? The glht function from multcomp can tell us. GLHT stands for General Linear Hypothesis Test.
library(multcomp)
result <- glht(pvcfit1, K)
result
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## Carl_vs_Bob == 0 -1.244
summary(result)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Carl_vs_Bob == 0 -1.2437 0.4048 -3.072 0.00391 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(result)
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Quantile = 2.0244
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Carl_vs_Bob == 0 -1.2437 -2.0632 -0.4243
glht can test multiple hypotheses at once. By default it applies a multiple testing correction when doing so. This is a generalization of Tukey’s Honestly Significant Differences.
K <- rbind(
Bob_vs_Alice = c(0, 1,0, 0,0,0,0,0,0,0),
Carl_vs_Alice = c(0, 0,1, 0,0,0,0,0,0,0),
Carl_vs_Bob = c(0, -1,1, 0,0,0,0,0,0,0))
result <- glht(pvcfit1, K)
summary(result)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Bob_vs_Alice == 0 -0.2625 0.4048 -0.648 0.79432
## Carl_vs_Alice == 0 -1.5063 0.4048 -3.721 0.00186 **
## Carl_vs_Bob == 0 -1.2437 0.4048 -3.072 0.01076 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(result)
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Quantile = 2.4378
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Bob_vs_Alice == 0 -0.2625 -1.2493 0.7243
## Carl_vs_Alice == 0 -1.5063 -2.4931 -0.5194
## Carl_vs_Bob == 0 -1.2437 -2.2306 -0.2569
plot(result)

We can also turn off the multiple testing correction.
summary(result, test=adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Bob_vs_Alice == 0 -0.2625 0.4048 -0.648 0.52059
## Carl_vs_Alice == 0 -1.5063 0.4048 -3.721 0.00064 ***
## Carl_vs_Bob == 0 -1.2437 0.4048 -3.072 0.00391 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
A reasonable compromise between these extremes is Benjamini and Hochberg’s False Discovery Rate (FDR) correction.
summary(result, test=adjusted("fdr"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = psize ~ operator + resin, data = pvc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Bob_vs_Alice == 0 -0.2625 0.4048 -0.648 0.52059
## Carl_vs_Alice == 0 -1.5063 0.4048 -3.721 0.00192 **
## Carl_vs_Bob == 0 -1.2437 0.4048 -3.072 0.00587 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
Finally, we can ask if any of the linear combinations is non-zero, i.e. whether the model with all three constraints applied can be rejected. This is equivalent to the anova( ) tests we have done earlier. (Note that while we have three constraints, the degrees of freedom reduction is 2, as any 2 of the constraints are sufficient. This makes me uneasy as it is reliant on numerical accuracy, better to just use any two of the constraints.)
summary(result, test=Ftest())
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## Bob_vs_Alice == 0 -0.2625
## Carl_vs_Alice == 0 -1.5063
## Carl_vs_Bob == 0 -1.2437
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 7.902 2 38 0.00135
pvcfit0 <- lm(psize ~ resin, data=pvc)
anova(pvcfit0, pvcfit1)
## Analysis of Variance Table
##
## Model 1: psize ~ resin
## Model 2: psize ~ operator + resin
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 40 70.533
## 2 38 49.815 2 20.718 7.902 0.00135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This demonstrates that the two methods of testing hypotheses–with the ANOVA test and with linear hypotheses–are equivalent.
An easier way to specify contrasts
So far we have been manually constructing linear hypotheses. The multcomp package automates this for some common situations. To compare all pairs of levels within a factor, use the mcp function, giving the factor to test as the name of an argument and specifying to test “Tukey” contrasts:
result <- glht(pvcfit1, mcp(resin="Tukey"))
To compare the first level of a factor to all other levels, specify to test “Dunnett” contrasts:
result <- glht(pvcfit1, mcp(resin="Dunnett"))
The linear hypotheses actually used can be inspected with result$linfct.
The emmeans package also automates many common comparisons. In particular, if you are working with models including interactions this package provides results for an “average” individual when examining main effects. By default it treats each level of other factors as being equally likely when calculating this average. It’s up to you to decide if this is sensible!
library(emmeans)
emmeans(pvcfit1, ~ resin)
emmeans(pvcfit1, pairwise ~ resin)
emmeans(pvcfit1, trt.vs.ctrl ~ resin)